Precision Forecasts of Land Cover Change

Chesapeake Watershed

Return to MUSA 801 Projects Page

This project was completed for the MUSA/Smart Cities Practicum course (MUSA 801) instructed by Michael Fichman and Matthew Harris. We are grateful to our instructors for their continued support and feedback. We would like to give special thanks to KC Filippino and Ben McFarlane from Hampton Roads Planning District Commission, and Dexter Locke from the United States Forest Service for providing data, insight, and support throughout the semester. This project would not have been possible without them.

View Dashboard

1.Introduction

1.1 Abstract

This project aims to develop a precision forecast model for land cover change at the Chesapeake Watershed. By leveraging high-resolution longitudinal land cover data provided by the Chesapeake Conservancy, the model will predict land cover conversions from pervious to impervious surfaces. This forecast will enable land use and environmental planners to identify where urban growth will occur, propose green infrastructure accordingly, and prioritize lands for protection. The model will be generalizable to the county level, incorporating only widely available inputs, thus allowing any municipality within the Chesapeake basin to replicate the analysis. This proof-of-concept project will demonstrate the utility of precision conservation in land protection and green infrastructure planning and provide a valuable tool for planners and policymakers across the region.

1.2 Background

The Chesapeake Bay watershed is an ecologically and economically significant resource, encompassing diverse ecosystems and supporting a multitude of industries, including agriculture, tourism, and fisheries. However, the region is facing increasing environmental challenges due to the combined effects of sea-level rise and land subsidence. As a result, the area has become the second-most vulnerable region in the nation to flooding and storm surge, only after New Orleans. Predicting land cover changes, particularly the conversion from pervious to impervious surfaces, is crucial in addressing these challenges and informing climate adaptation and mitigation planning.

Our project focuses on three distinct counties within the Chesapeake Bay watershed, representing varying development contexts. Portsmouth is the urban prototype characterized by its dense residential, commercial, and industrial areas. James City County exemplifies a suburban context, with a mix of rural, suburban, and urban development and a diverse landscape encompassing forests, wetlands, and historic sites. Lastly, Isle of Wight County represents the rural aspect, predominantly characterized by agriculture, forestry, and extensive natural habitats. By considering these diverse counties, we can develop a comprehensive and generalizable model to predict land cover changes across various regional development scenarios.

1.3 Motivation & Use Case

Building resilient communities is a top priority for the Hampton Roads Planning District Commission (HRPDC). To support this goal, the HRPDC has established a green infrastructure plan that seeks to identify and prioritize a network of valuable conservation lands. This plan aims to achieve multiple benefits, such as habitat protection, drinking water supply protection, stormwater management, and recreational opportunities.

A crucial component of this plan involves developing a model to forecast potential future growth and identify areas of the green infrastructure network that are most at risk for development. Our project aims to create a forecast that enables land use and environmental planners to pinpoint where urban growth is likely to occur, propose green infrastructure accordingly, and prioritize lands for protection.

This proof-of-concept project demonstrates the utility of precision conservation in land protection and green infrastructure planning, providing a valuable tool for planners and policymakers across the region. For example, Andrew, the head of the Green Infrastructure Team from Chesapeake Conservancy, and his team will use our web app to make informed decisions on which regions have the highest priority to receive funding.

2. Data and Methods

2.1 Understanding Land Cover data

The Chesapeake Conservancy provided high-resolution land cover data for 2013/14 and 2017/18. This vast raster dataset boasts an impressive 1-meter accuracy, offering 900 times more detail than the commonly used 30-meter resolution National Land Cover Dataset. Such a level of detail is crucial for capturing subtle changes in land cover.

The land cover classification includes pervious surfaces such as tree canopies and shrubs, which allow water to infiltrate the ground. In contrast, impervious surfaces encompass categories like roads and structures that prevent water infiltration, leading to increased runoff and potential flooding issues. Even though water and wetland are often considered impervious surfaces, in this study, we classify them as pervious surfaces due to their dynamic nature, interaction with groundwater, floodplain connectivity, and the critical functions of wetlands in water storage and infiltration.

Data Source

2.2 Ohter data

To understand how landcover change is affected by other environmental, social and economical factors, we also acquire data from the following source:

• DEM ( 1 arc-second) from USGS

• Soil Data from the Web Soil Survey

• census tract-level data (2014, 2018, 2021) from four years of the American Community Survey (ACS)

2.3 Unit of analysis

To facilitate our future analysis and ensure consistency, we gathered and resampled all data to a 10 x 10-meter resolution as the basic analysis unit. We opted for a 10-meter resolution as it provides sufficient detail for planning purposes while minimizing noise and reducing the dataset size, thereby enabling a faster calculation process.

3. Exploratory Analysis

3.1 Land Cover Type and Change Dynamics

We began by comparing the land cover patterns and changes across the three counties. This comparison helped us identify each county’s unique characteristics and trends, providing insights into how urban, suburban, and rural contexts affected land cover conversion. Land cover changes are highly related to the existing land cover type. For James City and Isle of Wight, tree canopy is the land cover type that undergoes the most change, while for Portsmouth, low vegetation experiences the most significant change.

3.2 Environmental Factors

Environmental factors played a significant role in land cover changes. We examined site characteristics that remained constant over time, such as elevation, slope, and soil type, to understand their influence on land cover conversion patterns. Understanding these static factors helped us predict how land cover might evolve in different environmental settings.

Take Portsmouth county as example, the plot below demonstrates that the mean slope in areas where land cover changed from pervious to impervious between 2014 and 2018 (1) is significantly higher than in areas that did not change or changed in the opposite direction (0).

3.3 Social and Economical Factors

Social and economic factors, which change over time, also impact land cover changes. We investigated variables such as population growth, economic development, and demographic changes to understand how these factors contributed to land cover conversion. By incorporating these dynamic factors into our analysis, we were able to better forecast land cover changes based on potential future scenarios.

The plots presented here illustrate that the percentage of white population change exhibits a different pattern in areas of land cover change compared to the overall county. This observation indicates that in Portsmouth County, land cover change is also associated with changes in census data, particularly the percentage of white population change.

4. Feature Engineering

4.1 Spatial Effects

Incorporating spatial relationships into our model was crucial for capturing the spatial influence of original land cover types on land cover change. We performed focal raster calculations and spatial lag calculations to create features that accounted for the spatial context of land cover conversion. These engineered features helped us develop a more accurate and robust predictive model, capable of capturing the nuances of land cover change across diverse development contexts.

4.2 Pairwise Correlations

To select variables for our predictive model, we visualize correlations between our numeric risk factors. This helps us ensure that our variables are not correlated with one another.

4.3 Final Features

Based on this plot, we select the following variables to use in our model: #### Dependent Variable * Whether the land cover change from impervious to previous

Independent Variables

• Existing land cover types
  • shrub
  • water
  • canopy
  • road
  • other (other impervious surface besides road)
  • impervious density
• Environmental factors
  • terrain (Elevation)
  • slope
  • soil type
• Socio-economic factors
  • population change
  • percent of white change
  • unit change
  • median household income change
  • area of census block group
• Spatial Lag factors
  • nonlinear transformation to predictor x
  • nonlinear transformation to predictor y

5. Modeling & Evaluation (Need Check!)

5.1 Model Building

To ensure computational efficiency and scalability, we downsampled the original dataset to 500,000 data points for model building and then fitted the selected model back to the whole dataset for future predictions. We employed geo cross-validation at the block group level to ensure the robustness of our model and avoid overfitting.

Performance evaluation was conducted using the confusion matrix for binary threshold setting and model selection. This approach allowed us to assess the accuracy, sensitivity, and specificity of the models and compare their performance.

5.2 Refine for better prediction

5.3 Model Type Selection

Recognizing that different types of models perform well on different datasets, we experimented with three model types for predicting land cover change: Random Forest, XGBoost, and Binomial Generalized Linear Model (GLM). These models were chosen due to their ability to handle complex interactions and non-linear relationships within the data.

After evaluating the performance of each model type, we selected Random Forest as the most suitable model for our analysis, as it demonstrated the best accuracy. Random Forest is an ensemble learning method that builds multiple decision trees and combines their results to improve overall accuracy and stability. This model is particularly well-suited for handling high-dimensional and noisy data, making it an ideal choice for predicting land cover changes in our study area.

5.4 Model Evaluation & Validation

roc

5.5 Comparison between Counties

The plots below indicate that spatial lag factors play a crucial role in land cover change prediction for all three counties. In addition, terrain, water, soil, and canopy emerge as important features across James City, Isle of Wight, and Portsmouth counties. For James City and Isle of Wight, slope is also a significant feature, whereas its importance is reduced in Portsmouth. This difference can be attributed to the urban nature of Portsmouth, where landforms less influence urban growth.

In Portsmouth, social and economic factors such as the percentage of white population change and median household income demonstrate greater significance in predicting land cover change. This suggests that urban growth in Portsmouth is more closely tied to socioeconomic factors, reflecting the distinctive characteristics of each county and the need for tailored land cover change modeling approaches.

6. Prediction & Error Analysis

In this section, we present the visualizations of our model’s results, highlighting areas where land cover changes from pervious to impervious surfaces are most likely to occur.

In James City County, the model predicts land cover changes are more likely to happen in areas with steeper slopes and in close proximity to existing development, particularly in the northern and southeastern parts of the county.

In Portsmouth County, the conversion from pervious to impervious surfaces is predicted to be more likely in two waterfront areas located in the eastern part of the county and in the midtown region.

In Isle of Wight County, land cover changes from pervious to impervious surfaces appear to be more scattered throughout the region, reflecting a more diverse pattern of potential development.

By comparing the predicted 2018 land cover change with the observed change, we can tell that errors in the prediction are likely to be around the true value, indicating the effectiveness of the model. Looking at the future prediction of 2021, we can see that many of the predicted changes are consistent with the observed change. However, the model has also identified new areas of land cover change. This suggests that the model is able to detect changes that were not previously observed, which can be valuable for land management and conservation efforts.

7. Conclusion

In conclusion, we successfully designed and set up an efficient workflow for testing and modeling land cover change predictions. Our chosen model, Random Forest, achieved a high balanced accuracy of approximately 98%, with performance metrics indicating that it generalizes well across different counties and land cover types. However, the model scripts require a considerable amount of time to run, which might pose challenges for replication and scalability.

Additionally, the model demonstrates lower accuracy in predicting land cover change compared to no change, which can be attributed to the uneven distribution of data. Despite these limitations, we are confident that our model and accompanying app can significantly contribute to the utility of precision conservation in land protection and green infrastructure planning. Ultimately, our project provides a valuable tool for planners and policymakers across the region, enabling them to make informed decisions on land use management and conservation strategies.

8. Code Appendix

########### Required Packages ###########
packages = c("stars","starsExtra", "abind","tigris","tidycensus","dyplyr",
             "bayesplot", "lme4","RcppEigen","knitr","yardstick",
             "tidyverse", "tidyr", "broom", "caret", "dials", "doParallel", "e1071", "earth",
             "ggrepel", "glmnet", "ipred", "klaR", "kknn", "pROC", "rpart", "randomForest",
             "sessioninfo", "tidymodels","ranger", "recipes", "workflows", "themis","xgboost",
             "sf", "nngeo", "mapview","raster")


########### Define Functions ###########

Load.fun <- function(x) { 
  x <- as.character(x) 
  if(isTRUE(x %in% .packages(all.available=TRUE))) { 
    eval(parse(text=paste("require(", x, ")", sep=""))) 
    print(paste(c(x, " : already installed; requiring"), collapse=''))
  } else { 
    #update.packages()
    print(paste(c(x, " : not installed; installing"), collapse=''))
    eval(parse(text=paste("install.packages('", x, "')", sep=""))) 
    print(paste(c(x, " : installed and requiring"), collapse=''))
    eval(parse(text=paste("require(", x, ")", sep=""))) 
  } 
} 

for(i in seq_along(packages)){
  packge <- as.character(packages[i])
  Load.fun(packge)
}

plotTheme <- function(base_size = 12, title_size = 12) {
  theme(
   text = element_text( color = "black"),
   plot.title = element_text(size = title_size, colour = "black"), 
        plot.subtitle = element_text(face="italic"),
       plot.caption = element_text(hjust=0),
         axis.ticks = element_blank(),
         panel.background = element_blank(),
         panel.grid.minor = element_line("grey90", size = 0.01),
         panel.grid.major = element_line("grey90", size = 0.01),
         panel.border = element_rect(colour = "white" , fill=NA, size=0),
         strip.background = element_rect(color = "white",fill='white'),
         strip.text = element_text(size=12),
         axis.title = element_text(size=12),
         axis.text = element_text(size=10),
         plot.background = element_blank(),
         legend.background = element_blank(),
         legend.title = element_text(colour = "black", face = "italic"),
         legend.text = element_text(colour = "black", face = "italic"),
         strip.text.x = element_text(size = 8)
       )
   }
mapTheme <- theme(plot.title =element_text(size=12),
                  plot.subtitle = element_text(size=8),
                  plot.caption = element_text(size = 6),
                  axis.line=element_blank(),
                  axis.text.x=element_blank(),
                  axis.text.y=element_blank(),
                  axis.ticks=element_blank(),
                  axis.title.x=element_blank(),
                  axis.title.y=element_blank(),
                  panel.background=element_blank(),
                  panel.border=element_blank(),
                  panel.grid.major=element_line(colour = 'transparent'),
                  panel.grid.minor=element_blank(),
                  legend.direction = "vertical", 
                  legend.position = "right",
                  plot.margin = margin(1, 1, 1, 1, 'cm'),
                  legend.key.height = unit(1, "cm"), legend.key.width = unit(0.2, "cm"))

palette2 <- c("#41b6c4","#253494")
palette4 <- c("#a1dab4","#41b6c4","#2c7fb8","#253494")
palette5 <- c("#ffffcc","#a1dab4","#41b6c4","#2c7fb8","#253494")
palette10 <- c("#f7fcf0","#e0f3db","#ccebc5","#a8ddb5","#7bccc4",
                        "#4eb3d3","#2b8cbe","#0868ac","#084081","#f7fcf0")

#########################################################################################################################################
###################################################  --------------------  ##############################################################
################################################### |                    | ##############################################################
################################################### |  DATA PREPARATION  | ##############################################################
################################################### |                    | ##############################################################
###################################################  --------------------  ##############################################################
#########################################################################################################################################

# load all the land cover data
port_14 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/lc/port_51740_lc_2014/port_51740_landcover_2014.tif") 
port_18 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/lc/port_51740_lc_2018/port_51740_landcover_2018.tif")
james_14 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/lc/jame_51095_lc_2014/jame_51095_landcover_2014.tif")
james_18 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/lc/jame_51095_lc_2018/jame_51095_landcover_2018.tif")

isle_14_10 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/lc/islelc_14_10x10.tif")
isle_18_10 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/lc/islelc_18_10x10.tif")

# load all the boundary data 
port_area <- st_read("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/boundary/portsmouth.shp")
james_area <- st_read("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/boundary/James.shp")
isle_area <- st_read("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/boundary/Isle_of_Wight.shp")

#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# | LANDCOVER DATASET  | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################


# Resample the raster image to 10x10 cells
port_14_10 <- st_warp(port_14, cellsize = 10, crs = st_crs(port_14))
port_18_10 <- st_warp(port_18, cellsize = 10, crs = st_crs(port_18))
jame_14_10 <- st_warp(jame_14, cellsize = 10, crs = st_crs(jame_14))
jame_18_10 <- st_warp(jame_18, cellsize = 10, crs = st_crs(jame_18))

port_area <- port_area %>% 
  st_transform(crs = st_crs(port_14))
jame_area <- jame_area %>% 
  st_transform(crs = st_crs(jame_14))
isle_area <- isle_area %>% 
  st_transform(crs = st_crs(isle_14))

# Crop the resampled raster image to the test area
port_14_10_crop <- st_crop(port_14_10, port_area)
port_18_10_crop <- st_crop(port_18_10, port_area)
jame_14_10_crop <- st_crop(jame_14_10, jame_area)
jame_18_10_crop <- st_crop(jame_18_10, jame_area)
isle_14_10_crop <- st_crop(isle_14, isle_area)
isle_18_10_crop <- st_crop(isle_18, isle_area)


# Create land cover change feature
port_change <- c(port_18_10_crop - port_14_10_crop)%>%
  mutate(lcchange = case_when(port_51740_landcover_2018.tif != 0 ~ 1,
                              port_51740_landcover_2018.tif != 0 ~ 0))
jame_change <- c(jame_18_10_crop - jame_14_10_crop)%>%
  mutate(lcchange = case_when(jame_51095_landcover_2018.tif != 0 ~ 1,
                              jame_51095_landcover_2014.tif!= 0 ~ 0))
isle_change <- c(isle_18_10_crop - isle_14_10_crop)%>%
  mutate(lcchange = case_when(islelc_18_10x10.tif != 0 ~ 1,
                              islelc_18_10x10.tif != 0 ~ 0))

# Reclassify the cropped raster image
port_14_10_rc <- port_14_10_crop %>% 
  mutate(
    originallc =  port_51740_landcover_2014.tif,
    lc = case_when(
    port_51740_landcover_2014.tif < 6  ~ 0, # previous
    port_51740_landcover_2014.tif >= 6 ~ 1 # impervious
  ))
port_18_10_rc <- port_18_10_crop %>% 
  mutate(
    originallc =  port_51740_landcover_2018.tif,
    lc = case_when(
    port_51740_landcover_2018.tif < 6  ~ 0,
    port_51740_landcover_2018.tif >= 6 ~ 1
  ))
jame_14_10_rc <- jame_14_10_crop %>% 
  mutate(originallc = jame_51095_landcover_2014.tif,
         lc = case_when(
         jame_51095_landcover_2014.tif < 6  ~ 0, # permeable
         jame_51095_landcover_2014.tif >= 6 ~ 1 # impermeable
         ))
jame_18_10_rc <- jame_18_10_crop %>% 
  mutate(originallc = jame_51095_landcover_2018.tif,
         lc = case_when(
           jame_51095_landcover_2018.tif < 6  ~ 0,
           jame_51095_landcover_2018.tif >= 6 ~ 1
         ))

isle_14_10_rc <- isle_14_10_crop %>% 
  mutate(originallc = case_when(
    islelc_14_10x10.tif == 1  ~ 1,
    islelc_14_10x10.tif == 2 ~ 2,
    islelc_14_10x10.tif == 3 ~ 3,
    islelc_14_10x10.tif == 4 ~ 4,
    islelc_14_10x10.tif == 5 ~ 5,
    islelc_14_10x10.tif == 6 ~ 6,
    islelc_14_10x10.tif == 7 ~ 7,
    islelc_14_10x10.tif == 8 ~ 8,
    islelc_14_10x10.tif == 9 ~ 9,
    islelc_14_10x10.tif == 10 ~ 10,
    islelc_14_10x10.tif == 11  ~ 11,
    islelc_14_10x10.tif == 12 ~ 12
  ),
  lc = case_when(
    islelc_14_10x10.tif < 6  ~ 0, # permeable
    islelc_14_10x10.tif >= 6 ~ 1 # impermeable
  ))
isle_18_10_rc <- isle_18_10_crop %>% 
  mutate(
    originallc =case_when(
      islelc_18_10x10.tif == 1  ~ 1,
      islelc_18_10x10.tif == 2 ~ 2,
      islelc_18_10x10.tif == 3 ~ 3,
      islelc_18_10x10.tif == 4 ~ 4,
      islelc_18_10x10.tif == 5 ~ 5,
      islelc_18_10x10.tif == 6 ~ 6,
      islelc_18_10x10.tif == 7 ~ 7,
      islelc_18_10x10.tif == 8 ~ 8,
      islelc_18_10x10.tif == 9 ~ 9,
      islelc_18_10x10.tif == 10 ~ 10,
      islelc_18_10x10.tif == 11  ~ 11,
      islelc_18_10x10.tif == 12 ~ 12
    ),
    lc = case_when(
      islelc_18_10x10.tif < 6  ~ 0,
      islelc_18_10x10.tif >= 6 ~ 1
    ))
# Create reclassed land cover change feature
port_change_rc <- (port_18_10_rc - port_14_10_rc)
jame_change_rc <- (jame_18_10_rc - jame_14_10_rc)
isle_change_rc <- (isle_18_10_rc - isle_14_10_rc)

#remove
rm(port_14_10,port_18_10,port_14,port_18)
rm(jame_14_10,jame_18_10,jame_14,jame_18)
rm(isle_14_10,isle_18_10,isle_14,isle_18)



#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |  CENSCUS DATASET   | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################


## load census

options(tigris_use_cache = TRUE)

census_api_key("bd4ce20561125a8480632db6acb29869e040ed08",install = TRUE)

# Define function to retrieve census data for a given county and year
get_county_data <- function(state, county, year) {
  get_acs(
    geography = "block group",
    variables = c("B01003_001E","B02001_002E","B19013_001E","B25002_001E","B06012_002E","B27011_008E"),
    year = year,
    state = state,
    county = county,
    geometry = TRUE,
    output = "wide"
  ) %>%
    st_transform(st_crs(port_18)) %>%
    rename(
      TotalPop = B01003_001E,
      Whites = B02001_002E,
      MedHHInc = B19013_001E,
      TotalUnit = B25002_001E
    ) %>%
    mutate(
      pctWhite = ifelse(TotalPop > 0, Whites / TotalPop * 100, 0),
      area = st_area(geometry)
    ) %>%
    dplyr::select(-NAME, -starts_with("B"),-Whites)%>%
    fill(everything())
}

# Define function to calculate changes in census variables between two years
calculate_changes <- function(data1, data2) {
  data2 %>%
    st_join(data1,left = TRUE) %>%
    mutate(
      popchange = (TotalPop.y - TotalPop.x) / (TotalPop.x * area.x),
      pctwhitechange = (pctWhite.y - pctWhite.x),
      Unitchange = (TotalUnit.y - TotalUnit.x) / area.x,
      MedHHIncchange = (MedHHInc.y - MedHHInc.x) / area.x,
      Area = area.x,
      GEOID = GEOID.x
    )%>%
    dplyr::select(-ends_with(".x"),-ends_with(".y"))
}

# Retrieve census data for Portsmouth in 2014 and 
portsmouth_14 <- get_county_data(state = "51", county = "Portsmouth", year = 2014)
portsmouth_18 <- get_county_data(state = "51", county = "Portsmouth", year = 2018)
portsmouth_21 <- get_county_data(state = "51", county = "Portsmouth", year = 2021)
# Calculate changes in census variables between 2014 and 2018
port_tract <- calculate_changes(data1 = portsmouth_14, data2 = portsmouth_18)
port_tract2 <- calculate_changes(data1 = portsmouth_18, data2 = portsmouth_21)

# Repeat for James City and Isle of Wight counties
james_city_14 <- get_county_data(state = "51", county = "James City", year = 2014)
james_city_18 <- get_county_data(state = "51", county = "James City", year = 2018)
james_city_21 <- get_county_data(state = "51", county = "James City", year = 2021)
jame_tract <- calculate_changes(data1 = james_city_14, data2 = james_city_18)
jame_tract2 <- calculate_changes(data1 = james_city_18, data2 = james_city_21)

isle_of_wight_14 <- get_county_data(state = "51", county = "Isle of Wight", year = 2014)
isle_of_wight_18 <- get_county_data(state = "51", county = "Isle of Wight", year = 2018)
isle_of_wight_21 <- get_county_data(state = "51", county = "Isle of Wight", year = 2021)
isle_tract <- calculate_changes(data1 = isle_of_wight_14, data2 = isle_of_wight_18)
isle_tract2 <- calculate_changes(data1 = isle_of_wight_18, data2 = isle_of_wight_21)

#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |    DEM DATASET     | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################


dem1 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/DEM/USGS_1_n38w077_20170509.tif")
dem2 <- read_stars("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/DEM/USGS_1_n37w077_20160315.tif")
dem <- st_mosaic(dem1, dem2)
dem2 <- st_warp(dem2, crs=st_crs(port_14))
dem_port <- st_crop (dem2, port_area)
dem_port <- st_warp(dem_port, port_14)
port_slope<- slope(dem_port)

dem1 <- st_warp(dem1, crs=st_crs(jame_area))
dem_james <- st_crop (dem1, jame_area)
dem_james <- st_warp(dem_james, jame_18_10_crop)
james_slope<- slope(dem_james)

dem2 <- st_warp(dem, crs=st_crs(isle_14_10_rc))
dem_isle <- st_crop (dem2, isle_area)
dem_isle <- st_warp(dem_isle, isle_14_10_crop)
isle_slope<- slope(dem_isle)
rm(dem1,dem2,dem)

#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |    SOIL DATASET    | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################

# Portmouth
soil <- st_read('~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/Soil/Port/port_soil.shp', crs= 'EPSG:4326')
soilclass <- soil%>%
  mutate(soil = case_when(
    Port__Rati == 'D' ~ 6,
    Port__Rati == 'C' ~ 5,
    Port__Rati == 'B/D' ~ 4,
    Port__Rati == 'B' ~ 3,
    Port__Rati == 'A/D' ~ 2,
    Port__Rati == 'A' ~ 1,
    is.na(Port__Rati) ~ 0
  )) %>%
  dplyr::select(-colnames(soil),geometry)%>%
  st_transform(st_crs(port_14_10_rc))

soil_port <- st_crop(soilclass, port_area)

# James City
soil <- st_read('~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/Soil/James/James_soil.shp', crs= 'EPSG:4326')
soilclass <- soil%>%
  mutate(soil = case_when(
    Report___8 == 'D' ~ 6,
    Report___8 == 'C' ~ 5,
    Report___8 == 'B/D' ~ 4,
    Report___8 == 'B' ~ 3,
    Report___8 == 'A/D' ~ 2,
    Report___8 == 'A' ~ 1,
    is.na(Report___8) ~ 0
  )) %>%
  dplyr::select(-colnames(soil),geometry)%>%
  st_transform(st_crs(jame_14_10_rc))

soil_jame <- st_crop(soilclass, jame_area)

# Isle of Wight
soil <- st_read('~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/Soil/isle/isle_soil.shp', crs= 'EPSG:4326')
soilclass <- soil%>%
  mutate(soil = case_when(
    isle__Ra_1 == 'D' ~ 6,
    isle__Ra_1 == 'C' ~ 5,
    isle__Ra_1 == 'B/D' ~ 4,
    isle__Ra_1 == 'B' ~ 3,
    isle__Ra_1 == 'A/D' ~ 2,
    isle__Ra_1 == 'A' ~ 1,
    is.na(isle__Ra_1) ~ 0
  )) %>%
  dplyr::select(-colnames(soil),geometry)%>%
  st_transform(st_crs(isle_14_10_rc))

soil_isle <- st_crop(soilclass, isle_area)
#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# | SPATIAL CALCULATION | #######################################################
#############################################  --------------------  ########################################################
#############################################################################################################################

port_14_pcnt_imperv = focal2(port_14_10_rc["lc"], matrix(1, 3, 3), "mean")
port_18_pcnt_imperv = focal2(port_18_10_rc["lc"], matrix(1, 3, 3), "mean")
jame_14_pcnt_imperv = focal2(jame_14_10_rc["lc"], matrix(1, 3, 3), "mean")
jame_18_pcnt_imperv = focal2(jame_18_10_rc["lc"], matrix(1, 3, 3), "mean")
isle_14_pcnt_imperv = focal2(isle_14_10_rc["lc"], matrix(1, 3, 3), "mean")
isle_18_pcnt_imperv = focal2(isle_18_10_rc["lc"], matrix(1, 3, 3), "mean")


# Function to calculate focal mean
calc_focal_mean <- function(df, field, ksize) {
  focal2(df[field], matrix(1, ksize, ksize), "mean") %>%
    focal2(., matrix(1, ksize, ksize), "mean") %>%
    focal2(., matrix(1, ksize, ksize), "mean")
}


#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |     PORTMOUTH      | #######################################################
#############################################  --------------------  ########################################################
#############################################################################################################################


# Road
port_14_road <- port_14_10_crop %>%
  mutate(road = as.numeric(port_51740_landcover_2014.tif == 9))

port_18_road <- port_18_10_crop %>%
  mutate(road = as.numeric(port_51740_landcover_2018.tif == 9))

road_14 <- calc_focal_mean(port_14_road, "road", 3)
road_18 <- calc_focal_mean(port_18_road, "road", 3)

# Canopy
port_14_canopy <- port_14_10_crop %>%
  mutate(canopy = as.numeric(port_51740_landcover_2014.tif == 3))

port_18_canopy <- port_18_10_crop %>%
  mutate(canopy = as.numeric(port_51740_landcover_2018.tif == 3))

canopy_14_mean <- calc_focal_mean(port_14_canopy, "canopy", 3)
canopy_18_mean <- calc_focal_mean(port_18_canopy, "canopy", 3)

# Shrub
port_14_shrub <- port_14_10_crop %>%
  mutate(shrub = as.numeric(port_51740_landcover_2014.tif == 2))

port_18_shrub <- port_18_10_crop %>%
  mutate(shrub = as.numeric(port_51740_landcover_2018.tif == 2))

shrub_14_mean <- calc_focal_mean(port_14_shrub, "shrub", 3)
shrub_18_mean <- calc_focal_mean(port_18_shrub, "shrub", 3)

# Water
port_14_water <- port_14_10_crop %>%
  mutate(water = as.numeric(port_51740_landcover_2014.tif == 1))

port_18_water <- port_18_10_crop %>%
  mutate(water = as.numeric(port_51740_landcover_2018.tif == 1))

water_14_mean <- calc_focal_mean(port_14_water, "water", 25)
water_18_mean <- calc_focal_mean(port_18_water, "water", 25)

# Other
port_14_other <- port_14_10_crop %>%
  mutate(other = as.numeric(port_51740_landcover_2014.tif == 8))

port_18_other <- port_18_10_crop %>%
  mutate(other = as.numeric(port_51740_landcover_2018.tif == 8))

other_14_mean <- focal2(port_14_other["other"], matrix(1, 3, 3), "mean")
other_14_mean <- focal2(other_14_mean["other"], matrix(1, 3, 3), "mean")
other_18_mean <- focal2(port_18_other["other"], matrix(1, 3, 3), "mean")
other_18_mean<- focal2(other_18_mean["other"], matrix(1, 3, 3), "mean")

## join all the data

port14 <- cbind(
  as.data.frame(canopy_14_mean)['canopy'],
  as.data.frame(road_14)['road'],
  as.data.frame(other_14_mean)['other'],
  as.data.frame(shrub_14_mean)['shrub'],
  as.data.frame(water_14_mean)['water'],
  as.data.frame(port_14_pcnt_imperv)['lc'] %>%
    rename(imperv = lc ),
  as.data.frame(port_change_rc)[4]%>%
    rename(lcchange = lc )
)


port18 <- cbind(
  as.data.frame(canopy_18_mean)['canopy'],
  as.data.frame(road_18)['road'],
  as.data.frame(other_18_mean)['other'],
  as.data.frame(shrub_18_mean)['shrub'],
  as.data.frame(water_18_mean)['water'],
  as.data.frame(port_18_pcnt_imperv)['lc'] %>%
    rename(imperv = lc )
)


port_14 <- 
  st_join(port_14_10_rc,port_tract) %>%
  st_join(.,dem_port)%>% # not changing
  st_join(.,port_slope) %>%
  st_join(., soil_port)


port_18 <- 
  st_join(port_18_10_rc,port_tract2) %>%
  st_join(.,dem_port)%>% # not changing
  st_join(.,port_slope) %>%
  st_join(., soil_port)

rm(dem_port,port_slope, port_14_10_crop, port_14_10_rc, port_14_canopy, port_14_other, port_14_road, port_14_water, port_18_10_crop, port_18_10_rc, port_18_canopy, port_18_other, port_18_road, port_18_water)


port14_df <- 
  as.data.frame(port_14) %>% 
  rename(
    terrain = USGS_1_n37w077_20160315.tif )%>%
  cbind(., port14)%>% na.omit()
port18_df <- 
  as.data.frame(port_18) %>% 
  rename(
    terrain = USGS_1_n37w077_20160315.tif) %>%
  cbind(., port18)%>%na.omit()

saveRDS(port14_df, "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/port14_df.rds")
saveRDS(port18_df, "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/port18_df.rds")


#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |      JAME CITY     | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################

# Road
jame_14_road <- jame_14_10_crop %>%
  mutate(road = as.numeric(jame_51095_landcover_2014.tif == 9))

jame_18_road <- jame_18_10_crop %>%
  mutate(road = as.numeric(jame_51095_landcover_2018.tif == 9))

road_14 <- calc_focal_mean(jame_14_road, "road", 3)
road_18 <- calc_focal_mean(jame_18_road, "road", 3)

# Canopy
jame_14_canopy <- jame_14_10_crop %>%
  mutate(canopy = as.numeric(jame_51095_landcover_2014.tif == 3))

jame_18_canopy <- jame_18_10_crop %>%
  mutate(canopy = as.numeric(jame_51095_landcover_2018.tif == 3))

canopy_14_mean <- calc_focal_mean(jame_14_canopy, "canopy", 3)
canopy_18_mean <- calc_focal_mean(jame_18_canopy, "canopy", 3)

# Shrub
jame_14_shrub <- jame_14_10_crop %>%
  mutate(shrub = as.numeric(jame_51095_landcover_2014.tif == 2))

jame_18_shrub <- jame_18_10_crop %>%
  mutate(shrub = as.numeric(jame_51095_landcover_2018.tif == 2))

shrub_14_mean <- calc_focal_mean(jame_14_shrub, "shrub", 3)
shrub_18_mean <- calc_focal_mean(jame_18_shrub, "shrub", 3)

# Water
jame_14_water <- jame_14_10_crop %>%
  mutate(water = as.numeric(jame_51095_landcover_2014.tif == 1))

jame_18_water <- jame_18_10_crop %>%
  mutate(water = as.numeric(jame_51095_landcover_2018.tif == 1))

water_14_mean <- calc_focal_mean(jame_14_water, "water", 25)
water_18_mean <- calc_focal_mean(jame_18_water, "water", 25)

# Other
jame_14_other <- jame_14_10_crop %>%
  mutate(other = as.numeric(jame_51095_landcover_2014.tif == 8))

jame_18_other <- jame_18_10_crop %>%
  mutate(other = as.numeric(jame_51095_landcover_2018.tif == 8))

other_14_mean <- focal2(jame_14_other["other"], matrix(1, 3, 3), "mean")
other_14_mean <- focal2(other_14_mean["other"], matrix(1, 3, 3), "mean")
other_18_mean <- focal2(jame_18_other["other"], matrix(1, 3, 3), "mean")
other_18_mean<- focal2(other_18_mean["other"], matrix(1, 3, 3), "mean")

## join all the data

jame14 <- cbind(
  as.data.frame(canopy_14_mean)['canopy'],
  as.data.frame(road_14)['road'],
  as.data.frame(other_14_mean)['other'],
  as.data.frame(shrub_14_mean)['shrub'],
  as.data.frame(water_14_mean)['water'],
  as.data.frame(jame_14_pcnt_imperv)['lc'] %>%
    rename(imperv = lc ),
  as.data.frame(jame_change_rc)['lc']%>%
    rename(lcchange = lc )
)

jame18 <- cbind(
  as.data.frame(canopy_18_mean)['canopy'],
  as.data.frame(road_18)['road'],
  as.data.frame(other_18_mean)['other'],
  as.data.frame(shrub_18_mean)['shrub'],
  as.data.frame(water_18_mean)['water'],
  as.data.frame(jame_18_pcnt_imperv)['lc'] %>%
    rename(imperv = lc )
)
jame_14 <- 
  st_join(jame_14_10_rc,jame_tract) %>%
  st_join(.,dem_james)%>% # not changing
  st_join(.,james_slope) %>%
  st_join(., soil_jame)


jame_18 <- 
  st_join(jame_18_10_rc,jame_tract2) %>%
  st_join(.,dem_james)%>% # not changing
  st_join(.,james_slope) %>%
  st_join(., soil_jame)

rm(dem_jame,jame_slope, jame_14_10_crop, jame_14_10_rc, jame_14_canopy, jame_14_other, jame_14_road, jame_14_water, jame_18_10_crop, jame_18_10_rc, jame_18_canopy, jame_18_other, jame_18_road, jame_18_water)

jame14_df <- 
  as.data.frame(jame_14) %>% 
  rename(
    terrain = USGS_1_n38w077_20170509.tif )%>%
  cbind(., jame14)  %>%
  na.omit() 
jame18_df <- 
  as.data.frame(jame_18) %>% 
  rename(
    terrain = USGS_1_n38w077_20170509.tif) %>%
  cbind(., jame18)%>%
  na.omit()

saveRDS(jame14_df, "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/jame14_df.rds")
saveRDS(jame18_df, "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/jame18_df.rds")

#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |   ISLE of WIGHT    | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################

# Road
isle_14_road <- isle_14_10_crop %>%
  mutate(road = as.numeric(islelc_14_10x10.tif == 9))

isle_18_road <- isle_18_10_crop %>%
  mutate(road = as.numeric(islelc_18_10x10.tif == 9))

road_14 <- calc_focal_mean(isle_14_road, "road", 3)
road_18 <- calc_focal_mean(isle_18_road, "road", 3)

# Canopy
isle_14_canopy <- isle_14_10_crop %>%
  mutate(canopy = as.numeric(islelc_14_10x10.tif == 3))

isle_18_canopy <- isle_18_10_crop %>%
  mutate(canopy = as.numeric(islelc_18_10x10.tif == 3))

canopy_14_mean <- calc_focal_mean(isle_14_canopy, "canopy", 3)
canopy_18_mean <- calc_focal_mean(isle_18_canopy, "canopy", 3)

# Shrub
isle_14_shrub <- isle_14_10_crop %>%
  mutate(shrub = as.numeric(islelc_14_10x10.tif == 2))

isle_18_shrub <- isle_18_10_crop %>%
  mutate(shrub = as.numeric(islelc_18_10x10.tif == 2))

shrub_14_mean <- calc_focal_mean(isle_14_shrub, "shrub", 3)
shrub_18_mean <- calc_focal_mean(isle_18_shrub, "shrub", 3)

# Water
isle_14_water <- isle_14_10_crop %>%
  mutate(water = as.numeric(islelc_14_10x10.tif == 1))

isle_18_water <- isle_18_10_crop %>%
  mutate(water = as.numeric(islelc_18_10x10.tif == 1))

water_14_mean <- calc_focal_mean(isle_14_water, "water", 25)
water_18_mean <- calc_focal_mean(isle_18_water, "water", 25)

# Other
isle_14_other <- isle_14_10_crop %>%
  mutate(other = as.numeric(islelc_14_10x10.tif == 8))

isle_18_other <- isle_18_10_crop %>%
  mutate(other = as.numeric(islelc_18_10x10.tif == 8))

other_14_mean <- focal2(isle_14_other["other"], matrix(1, 3, 3), "mean")
other_14_mean <- focal2(other_14_mean["other"], matrix(1, 3, 3), "mean")
other_18_mean <- focal2(isle_18_other["other"], matrix(1, 3, 3), "mean")
other_18_mean<- focal2(other_18_mean["other"], matrix(1, 3, 3), "mean")

## join all the data

isle14 <- cbind(
  as.data.frame(canopy_14_mean)['canopy'],
  as.data.frame(road_14)['road'],
  as.data.frame(other_14_mean)['other'],
  as.data.frame(shrub_14_mean)['shrub'],
  as.data.frame(water_14_mean)['water'],
  as.data.frame(isle_14_pcnt_imperv)['lc'] %>%
    rename(imperv = lc ),
  as.data.frame(isle_change_rc)[c(5)]%>%
    rename(lcchange = lc )
)

isle18 <- cbind(
  as.data.frame(canopy_18_mean)['canopy'],
  as.data.frame(road_18)['road'],
  as.data.frame(other_18_mean)['other'],
  as.data.frame(shrub_18_mean)['shrub'],
  as.data.frame(water_18_mean)['water'],
  as.data.frame(isle_18_pcnt_imperv)['lc'] %>%
    rename(imperv = lc )
)

isle14 <-  cbind(isle_14_10_crop,canopy_14_mean, road_14, other_14_mean, shrub_14_mean, water_14_mean, isle_14_pcnt_imperv, isle_change_rc,isle_change) # change impervious(0, 1, -1) 1 is targeted
rm(canopy_14_mean, road_14, other_14_mean, shrub_14_mean, water_14_mean, isle_14_pcnt_imperv, isle_change_rc)

isle18 <- cbind(isle_18_10_crop, canopy_18_mean, road_18,other_18_mean, shrub_18_mean, water_18_mean,isle_18_pcnt_imperv)
rm(canopy_18_mean, road_18, other_18_mean, shrub_18_mean, water_18_mean, isle_18_pcnt_imperv)

isle_14 <- 
  st_join(isle_14_10_rc,isle_tract) %>%
  st_join(.,dem_isle)%>% # not changing
  st_join(.,isle_slope) %>%
  st_join(., soil_isle)


isle_18 <- 
  st_join(isle_18_10_rc,isle_tract2) %>%
  st_join(.,dem_isle)%>% # not changing
  st_join(.,isle_slope) %>%
  st_join(., soil_isle)

rm(dem_isle,isle_slope, isle_14_10_crop, isle_14_10_rc, isle_14_canopy, isle_14_other, isle_14_road, isle_14_shrub, isle_18_10_crop, isle_18_10_rc, isle_18_canopy, isle_18_other, isle_18_road, isle_18_shrub)

isle14_df <- 
  as.data.frame(isle_14)  %>%
  rename(terrain = USGS_1_n38w077_20170509.tif )%>%
  cbind(., isle14) %>% 
  na.omit()
isle18_df <- 
  as.data.frame(isle_18) %>% 
  rename(terrain = USGS_1_n38w077_20170509.tif) %>%
  cbind(., isle18)%>%
  na.omit()

saveRDS(isle14_df, "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/isle14_df.rds")
saveRDS(isle18_df, "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/isle18_df.rds")

#########################################################################################################################################
###################################################  --------------------  ##############################################################
################################################### |                    | ##############################################################
################################################### |      MODELING      | ##############################################################
################################################### |                    | ##############################################################
###################################################  --------------------  ##############################################################
#########################################################################################################################################



set.seed(717)


"%!in%" <- Negate("%in%")

isle14_df <- readRDS( "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/isle14_df.rds")
jame14_df <- readRDS( "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/jame14_df.rds")
port14_df <- readRDS( "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/port14_df.rds")



# Data preparation
data <- data %>% 
  mutate(
    lcre = case_when(lcchange == -1 ~ 0,
                     lcchange == 1 ~ 1,
                     lcchange == 0 ~ 0)
  ) %>% 
  mutate(
    lcchange = as.factor(lcchange),
    lcre = as.factor(lcre)
  )


# Initial split for training and test
data_split <- initial_split(data, strata = "lcchange", prop = 0.8)
data_train <- training(data_split)
data_test <- testing(data_split)

apply_model_rf <- function(data) {
  # Sample data
  balanced_data <- data_train[c(sample(which(data_train$lcre == 0), sum(data_train$lcre == 1) *10) , which(data_train$lcre == 1)),]
  # Cross-validation
  cv_splits_geo <- group_vfold_cv(balanced_data, group = "GEOID")
  
  # Create recipe
  model_rec <- recipe(lcre ~ ., data = balanced_data) %>%
    update_role(GEOID, new_role = "GEOID") %>% #78
    update_role(lcchange, new_role = "lcchange") %>%
    update_role(port_51740_landcover_2014.tif, new_role = "originallc") %>%
    update_role(lc, new_role = "lc")# %>%
  step_ns(x, y, options = list(df = 1))
  
  # Model specifications
  rf_plan <- rand_forest() %>%
    set_args(mtry  = tune()) %>%
    set_args(min_n = tune()) %>%  # 100,1000,seq=200
    set_args(trees =  tune()) %>% # 1000
    #set_args(max_depth  =  10) %>% 
    set_engine("ranger", importance = "impurity") %>% 
    set_mode("classification")
  
  rf_grid <- expand.grid(mtry = c(5,10,15),
                         min_n  = c(100,500,900),
                         trees =  c(1000,1500)
  )
  
  # Create the workflow
  rf_wf <-
    workflow() %>% 
    add_recipe(model_rec) %>% 
    add_model(rf_plan)
  
  # Tune hyperparameters
  control <- control_resamples(save_pred = TRUE, verbose = TRUE)
  metrics <- metric_set(f_meas)
  
  rf_tuned <- rf_wf %>%
    tune_grid(resamples = cv_splits_geo,
              grid      = rf_grid,
              control   = control,
              metrics   = metrics)
  
  # Select best model
  rf_best_params <- select_best(rf_tuned, metric = "f_meas")
  rf_best_wf <- finalize_workflow(rf_wf, rf_best_params)
  
  return(rf_best_wf)
}

# Fit multiple models
num_models <- 5
predict_dfs <- list()
rf_wfs <- list()
for (i in 1:num_models) {
  # Split data
  rf_wf <- apply_model_rf(data_train)
  
  # Make predictions
  rf_full_wf <- rf_wf %>%
    fit(data)
  predict_df <-
    predict(rf_full_wf , new_data = data, type = "prob")
  # Add to list
  predict_dfs[[i]] <- predict_df[2]
  rf_wfs[[i]] <- rf_full_wf
}


# Apply model to three different datasets
results_isle<- apply_model(isle14_df)
results_port<- apply_model(port14_df)
results_jame <- apply_model(jame14_df)

#########################################################################################################################################
###################################################  --------------------  ##############################################################
################################################### |                    | ##############################################################
################################################### |     PREDICTING     | ##############################################################
################################################### |                    | ##############################################################
###################################################  --------------------  ##############################################################
#########################################################################################################################################

isle18_df <- readRDS("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/isle18_df.rds")
port18_df <- readRDS("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/port18_df.rds")
jame18_df <- readRDS("~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/jame18_df.rds")

#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |      JAME CITY     | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################
# Fit the best model to the whole dataset
full_fit_rf <- results_jame[1][2]%>%
  fit(data = jame14_df)

# See the prediction for the full dataset
predict_df <- predict(full_fit_rf, new_data = jame14_df)
predict_df2 <- cbind(jame14_df,predict_df) %>%
  mutate(error = as.numeric(.pred_class) - as.numeric(lcre))

# Use the best model to predict for the future
predict_df3 <- predict(full_fit_rf, new_data = jame18_df %>% mutate(lcchange = as.factor(lc)))
predict_df3 <- cbind(jame18_df,predict_df3) 
saveRDS(predict_df3,'~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/predict_DF3_jame.rds')

# extract final model object
rf_full_mod <- extract_fit_parsnip(full_fit_rf)
saveRDS(rf_full_mod,'~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/best_rf1_jame.rds')

imp <- importance(rf_full_mod$fit, type = 1)
if (is.matrix(imp)) {
  imp_df <- data.frame(variable = rownames(imp), importance = imp[, 1])
} else {
  imp_df <- data.frame(variable = names(imp), importance = imp)
}


ggplot(imp_df, aes(x = reorder(variable, importance), y = importance, fill = "#41b6c4", alpha=0.5)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_identity() +
  labs(x = "Variable", y = "Importance") +
  ggtitle("Variable Importance")
#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |      PORTMOUTH     | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################

# Fit the best model to the whole dataset
full_fit_rf <- results_port[1][2]%>%
  fit(data = port14_df)

# See the prediction for the full dataset
predict_df <- predict(full_fit_rf, new_data = port14_df)
predict_df2 <- cbind(port14_df,predict_df) %>%
  mutate(error = as.numeric(.pred_class) - as.numeric(lcre))

# Use the best model to predict for the future
predict_df3 <- predict(full_fit_rf, new_data = port18_df %>% mutate(lcchange = as.factor(lc)))
predict_df3 <- cbind(port18_df,predict_df3) 
saveRDS(predict_df3,'~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/predict_DF3_port.rds')

# extract final model object
rf_full_mod <- extract_fit_parsnip(full_fit_rf)
saveRDS(rf_full_mod,'~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/best_rf1_port.rds')

#############################################################################################################################
#############################################  --------------------  ########################################################
############################################# |   ISLE OF WIGHT    | ########################################################
#############################################  --------------------  ########################################################
#############################################################################################################################
# Fit the best model to the whole dataset
full_fit_rf <- results_jame[1][2]%>%
  fit(data = port14_df)

# See the prediction for the full dataset
predict_df <- predict(full_fit_rf, new_data = isle14_df)
predict_df2 <- cbind(isle14_df,predict_df) %>%
  mutate(error = as.numeric(.pred_class) - as.numeric(lcre))

# Use the best model to predict for the future
predict_df3 <- predict(full_fit_rf, new_data = isle18_df %>% mutate(lcchange = as.factor(lc)))
predict_df3 <- cbind(isle18_df,predict_df3) 
saveRDS(predict_df3,'~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/predict_DF3_isle.rds')

# extract final model object
rf_full_mod <- extract_fit_parsnip(full_fit_rf)
saveRDS(rf_full_mod,'~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/best_rf1_isle.rds')

#########################################################################################################################################
###################################################  --------------------  ##############################################################
################################################### |                    | ##############################################################
################################################### |      MAPPING       | ##############################################################
################################################### |                    | ##############################################################
###################################################  --------------------  ##############################################################
#########################################################################################################################################


predict_df2<- readRDS('~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/predict_DF3_port.rds')
s <- st_as_stars(predict_df2,
                 dimensions=st_dimensions(
                   x=sort(unique(df$x)),
                   y=sort(unique(df$y)), 
                   point=TRUE),
                 dims=c('x','y'))
st_crs(s) <- st_crs(port_14)
s <- s%>%
  st_transform(4326)
x <- s['.pred_class']
x.or <-s['lcre']
error <- s['error']
thresh_value <- 1
x[x != thresh_value] <- NA
x_sf <- st_union(st_as_sf(x))
# write the sf object to a GeoJSON file
geojson_path <- "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/port_predict_xgb.geojson"
st_write(x_sf, geojson_path)



predict_df2<- readRDS('~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/predict_DF3_port.rds')
s <- st_as_stars(predict_df2,
                 dimensions=st_dimensions(
                   x=sort(unique(df$x)),
                   y=sort(unique(df$y)), 
                   point=TRUE),
                 dims=c('x','y'))
st_crs(s) <- st_crs(port_18)
s <- s%>%
  st_transform(4326)
x <- s['.pred_class']
x.or <-s['lcre']
error <- s['error']
thresh_value <- 1
x[x != thresh_value] <- NA
x_sf <- st_union(st_as_sf(x))
# write the sf object to a GeoJSON file
geojson_path <- "~/Github/Precision-Forecasts-of-Land-Cover-Change/Data/output/port_test_4326.geojson"
st_write(x_sf, geojson_path)